In [1]:
from netCDF4 import Dataset, num2date
from datetime import datetime
import glob
from mpl_toolkits.basemap import Basemap
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import calendar
from scipy.stats import describe
from scipy import interpolate
from tqdm import tqdm_notebook as tqdm
import oceansdb
#import seawater as sw
#from seawater.library import T90conv
In [22]:
def load_netcdf(path = "Argo_South_60"):
# load_netcdf function input: path - default is "Argo_South_60"
files = glob.glob("data/" + path + "/**/*.nc", recursive=True)
# glob.glob - returns a list of all the files in a given folder with the extension 'nc'
# recursive=true - keeps on looking for all the files till the end, or the last file that exist
# creating lists
data = {
"lat": [],
"lon": [],
"datetime": [],
"pressure": [],
"temperature": [],
"salinity": []
}
for f in tqdm(files):
# tqdm - make your loops show a progress meter
d = Dataset(f)
# DATASET - Creating/Opening/Closing a netCDF file.
#This is the method used to open an existing netCDF file.
lat = d.variables["LATITUDE"][:]
mask = lat < -60
lon = d.variables["LONGITUDE"][:]
data["lat"].extend(lat[mask])
data["lon"].extend(lon[mask])
#1 The extend () method adds the specified list elements to the end of the current list.
#2 This method does not return any value but add the content to existing list.
juld = d.variables["JULD"][:]
units = d.variables["JULD"].getncattr('units')
#1 getncattr - reads the attributes of given variable i.e units
dates = num2date(juld, units, "standard")
#1 num2date - converts the format of the loaded datetime to standard format YYYY-MM-dd
data["datetime"].extend(dates[mask])
data["pressure"].extend(d.variables["PRES_ADJUSTED"][:][mask])
data["temperature"].extend(d.variables["TEMP_ADJUSTED"][:][mask])
try:
data["salinity"].extend(d.variables["PSAL_ADJUSTED"][:][mask])
except:
data["salinity"].extend(np.full(len(mask), np.nan))
#1 array() - allows one to do operations on a list
#2 NumPy was imported as np
#3 If we want to do math on a homogeneous array of numeric data, then it is much better to use NumPy,
# which can automatically vectorize operations on complex multi-dimensional arrays
for k,v in data.items():
data[k] = np.array(v)
return data
def plot(lats, lons, z = [], title = "Argo profiles south of 60S", cbtitle = "Number of points in bin", vmax=None):
# setup north polar stereographic basemap.
# The longitude lon_0 is at 6-o'clock, and the
# latitude circle boundinglat is tangent to the edge
# of the map at lon_0. Default value of lat_ts
# (latitude of true scale) is pole.
fig = plt.figure(figsize=(15,15))
m = Basemap(projection='spstere',boundinglat=-55,lon_0=180,resolution='l')
m.drawcoastlines()
m.fillcontinents(color='black',lake_color='aqua')
# draw parallels and meridians.
m.drawparallels(np.arange(-60, 0, 20))
m.drawmeridians(np.arange(-180, 181, 20), labels=[1,1,0,1])
#m.drawmapboundary(fill_color='aqua')
plt.title("{} {}".format(len(lats), title))
x, y = m(lons, lats)
if len(z) == 0:
hh, locx, locy = np.histogram2d(x, y, bins=100)
# Sort the points by density, so that the densest points are plotted last
z = np.array([hh[np.argmax(a<=locx[1:]),np.argmax(b<=locy[1:])] for a,b in zip(x,y)])
#print(describe(z))
idx = z.argsort()
x, y, z = x[idx], y[idx], z[idx]
#m.imshow(heatmap, interpolation='bicubic', cmap="jet")
m.scatter(x, y, c=z, s=5, alpha=1, cmap="jet", vmax=vmax)
cb = plt.colorbar()
cb.set_label(cbtitle, rotation=270)
plt.show()
def plot_time(dts, title, label):
fig, ax = plt.subplots(1, 1, figsize=(20,10))
for i, dt in enumerate(dts):
#1 enumerate() method adds a counter to an iterable and returns it in a form of enumerate object.
# This enumerate object can then be used directly in for loops or be converted into a list of tuples using list() method.
k = label[i]
df = pd.DataFrame({k: dt})
#1 Pandas was imported as pd
#2 Dataframe class - class pandas.DataFrame(data=None, index=None, columns=None, dtype=None, copy=False)
# Two-dimensional size-mutable, potentially heterogeneous tabular data structure with labeled axes (rows and columns).
# Arithmetic operations align on both row and column labels.
# Can be thought of as a dict-like container for Series objects
df = df.groupby(by=df[k].dt.date).count()
#1 groupby - groups data by a certain specified variable type
df.plot(ax=ax, style='.', markersize=3)
ax.xaxis.set_major_locator(mdates.YearLocator())
ax.set_title(title, fontsize=18)
ax.set_xlabel("Time",fontsize=14)
ax.set_ylabel("Number of profiles",fontsize=14)
#ax.set_xlabel("")
ax.legend(fontsize=14)
plt.show()
In [ ]:
argo = load_netcdf()
In [30]:
plot(argo["lat"], argo["lon"], vmax=100)
sst = [x[0] for x in argo["temperature"]]
plot(argo["lat"], argo["lon"], sst, title = "Argo SST", cbtitle = "Temperature °C")
max_depth = [np.nanmax(x) for x in argo["pressure"]]
plot(argo["lat"], argo["lon"], max_depth, title = "Argo Max Depth", cbtitle = "Depth (decibar pressure)")
In [32]:
print(min(argo["datetime"]), max(argo["datetime"]))
plot_time([argo["datetime"]], "Argo float dailystamps", ["Argo"])
In [34]:
seal = load_netcdf("seal")
plot(seal["lat"], seal["lon"], title = "Seal data south of 60S", vmax=100)
sst = [x[0] for x in seal["temperature"]]
plot(seal["lat"], seal["lon"], sst, title = "Seal SST", cbtitle = "Temperature °C")
max_depth = [np.nanmax(x) for x in seal["pressure"]]
plot(seal["lat"], seal["lon"], max_depth, title = "Seal Max Depth", cbtitle = "Depth (decibar pressure)")
In [36]:
print(min(seal["datetime"]), max(seal["datetime"]))
plot_time([seal["datetime"]], "Seal data dailystamps", ["Elephant seal"])
In [38]:
all_lats = np.concatenate((seal["lat"], argo["lat"]))
all_lons = np.concatenate((seal["lon"], argo["lon"]))
plot(all_lats, all_lons, title="Argo + seal data south of 60S", vmax=100)
In [40]:
plot_time([seal["datetime"], argo["datetime"]], "Argo + seal data dailystamps", ["Seal", "Argo"])
In [165]:
def load_netcdf_grid(path = "Argo_South_60", resolution = 1, target_var = "temp", with_depth = False, filter_month = None, filter_season = None, filter_start = datetime(2000,1,1), filter_end = datetime(2020,1,1)):
if with_depth:
grid = np.full(shape=[15 * resolution, 360 * resolution, 200], fill_value=np.nan)
gridcounts = np.zeros(shape=[15 * resolution, 360 * resolution, 200])
else:
grid = np.full(shape=[15 * resolution, 360 * resolution], fill_value=np.nan)
files = glob.glob("data/{}/**/*.nc".format(path), recursive=True)
#1 glob.glob - returns a list of all the files in a given folder with the extension 'nc'
#2 recursive=true - keeps on looking for all the files till the end, or the last file that exist
for f in tqdm(files[:]):
# tqdm - make your loops show a progress meter
d = Dataset(f)
#1 DATASET - Creating/Opening/Closing a netCDF file.
# This is the method used to open an existing netCDF file.
lat = d.variables["LATITUDE"][:]
mask = (lat < -60) & (lat > -74.5)
juld = d.variables["JULD"][:]
units = d.variables["JULD"].getncattr('units')
dates = num2date(juld, units, "standard")
datemask = (dates > filter_start) & (dates < filter_end)
mask &= datemask
if filter_month:
datemask = np.array([d.month == filter_month for d in dates])
mask &= datemask
if filter_season:
if filter_season == "Summer":
datemask = np.array([d.month >= 10 or d.month <= 3 for d in dates])
elif filter_season == "Winter":
datemask = np.array([d.month > 3 and d.month < 10 for d in dates])
mask &= datemask
if any(mask):
lat = np.round((np.abs(lat[mask]) - 60) * resolution).astype(int)
lon = d.variables["LONGITUDE"][mask]
lon = np.round((lon + 180) * resolution).astype(int)
lon[lon == 360 * resolution] = 0
if with_depth:
pres = np.floor(d.variables["PRES_ADJUSTED"][mask] / 10).astype(int)
temp = d.variables["TEMP_ADJUSTED"][mask]
if target_var == "sal":
try:
sal = d.variables["PSAL_ADJUSTED"][mask]
except:
print("No salinity for {}".format(f))
continue
else:
pres = d.variables["PRES_ADJUSTED"][mask]
temp = d.variables["TEMP_ADJUSTED"][mask, 0]
for x in np.unique(lon):
for y in np.unique(lat):
if with_depth:
ptmask = (lon == x) & (lat == y)
for z in np.unique(pres[ptmask]):
if z>= 200 or np.ma.is_masked(z):
continue
depthmask = pres[ptmask] == z
if target_var == "temp":
values_at_pt = temp[ptmask][depthmask].compressed()
elif target_var == "sal":
values_at_pt = sal[ptmask][depthmask].compressed()
grid[y, x, z] = np.nansum((grid[y,x,z], np.nansum(values_at_pt)))
gridcounts[y, x, z] += len(values_at_pt)
else:
ptmask = (lon == x) & (lat == y)
if target_var == "n_point":
v = np.sum(ptmask)
if not np.ma.is_masked(v):
############# Nick see this line below added on 21/10/2019
if v != 0:
#1 added a if statement checking if number of points is greater than 0
grid[y, x] = np.nansum((grid[y, x], v))
elif target_var == "temp":
temps_at_pt = temp[ptmask]
v = np.ma.mean(temps_at_pt)
if not np.ma.is_masked(v):
grid[y, x] = np.nanmean((grid[y, x], v))
elif target_var == "depth":
pres_at_pt = pres[ptmask]
if len(pres_at_pt):
v = np.ma.max(pres[ptmask])
if not np.ma.is_masked(v):
grid[y, x] = np.nanmax((grid[y, x], v))
if with_depth:
grid /= gridcounts
print(np.nanmin(grid), np.nanmean(grid), np.nanmax(grid))
return grid
# This version seems a lot more inefficient for some reason
def load_netcdf_grid_new(platform = "argo", resolution = 1, target_var = "temp", filter_month = None, filter_season = None, filter_start = datetime(2000,1,1), filter_end = datetime(2020,1,1)):
grid = np.full(shape=[15 * resolution, 360 * resolution, 200], fill_value=np.nan)
if platform == "argo":
data = argo
else:
data = seal
mask = (data["lat"] < -60) & (data["lat"] > -74.5)
datemask = (data["datetime"] > filter_start) & (data["datetime"] < filter_end)
mask &= datemask
if filter_month:
datemask = np.array([d.month == filter_month for d in data["datetime"]])
mask &= datemask
if filter_season:
if filter_season == "Summer":
datemask = np.array([d.month >= 10 or d.month <= 3 for d in data["datetime"]])
elif filter_season == "Winter":
datemask = np.array([d.month > 3 and d.month < 10 for d in data["datetime"]])
mask &= datemask
if any(mask):
lat = np.round((np.abs(data["lat"][mask]) - 60) * resolution).astype(int)
lon = data["lon"][mask]
lon = np.round((lon + 180) * resolution).astype(int)
lon[lon == 360 * resolution] = 0
for y in tqdm(np.unique(lat)):
for x in tqdm(np.unique(lon)):
ptmask = (lon == x) & (lat == y)
depths = set()
for pres in data["pressure"][mask][ptmask]:
print(np.floor(pres / 10).compressed().astype(int))
depths.update(np.floor(pres / 10).compressed().astype(int))
print(len(depths))
for z in tqdm(depths):
depthmask = [np.floor(pres / 10) == z for pres in data["pressure"][mask][ptmask]]
values = []
for i in range(len(depthmask)):
values.extend(data["temperature"][mask][ptmask][i][depthmask[i]])
#print(values)
grid[y, x, z] = np.ma.mean(values)
return grid
def plot_grid(grid, resolution = 1, title="Argo data\n0-10 decibars mean temperature", cbtitle="Temperature °C", vmin=None, vmax=None):
fig = plt.figure(figsize=(15,15))
m = Basemap(projection='spstere',boundinglat=-55,lon_0=180,resolution='l')
m.drawcoastlines()
m.fillcontinents(color='black',lake_color='aqua')
# draw parallels and meridians.
m.drawparallels(np.arange(-60, -80, -1))
m.drawmeridians(np.arange(-180, 181, 20), labels=[1,1,0,1])
#m.drawmapboundary(fill_color='aqua')
plt.title("\n{} gridded at {}° resolution".format(title, 1/resolution),fontsize = 16)
x = np.arange(-180, 180.1, 1/resolution)
y = np.arange(-60, -75.1, 1/-resolution)
x, y = np.meshgrid(x, y)
x, y = m(x, y)
m.pcolormesh(x, y, grid, cmap="jet", vmin=vmin, vmax=vmax)
cb = plt.colorbar()
cb.set_label(cbtitle, rotation=270, fontsize = 12)
plt.show()
def plot_grid_a(grid, resolution = 1, title="Argo data\n0-10 decibars mean temperature", cbtitle="Temperature °C", vmin=None, vmax=None):
fig = plt.figure(figsize=(15,15))
m = Basemap(projection='spstere',boundinglat=-55,lon_0=180,resolution='l')
m.drawcoastlines()
m.fillcontinents(color='black',lake_color='aqua')
# draw parallels and meridians.
m.drawparallels(np.arange(-60, 0, 20))
m.drawmeridians(np.arange(-180, 181, 20), labels=[1,1,0,1])
#m.drawmapboundary(fill_color='aqua')
plt.title("{}".format(title),fontsize = 16)
#plt.title("{} at {}° resolution".format(title, 1/resolution),fontsize = 16)
x = np.arange(-180, 180.1, 1/resolution)
y = np.arange(-60, -75.1, 1/-resolution)
x, y = np.meshgrid(x, y)
x, y = m(x, y)
m.pcolormesh(x, y, grid, cmap="jet", vmin=vmin, vmax=vmax)
cb = plt.colorbar()
cb.set_label(cbtitle, rotation=270, fontsize = 12)
plt.show()
def interp_nans(grid, method='linear'):
x = np.arange(0, grid.shape[1])
y = np.arange(0, grid.shape[0])
mask = np.isnan(grid)
xx, yy = np.meshgrid(x, y)
x1 = xx[~mask]
y1 = yy[~mask]
newarr = grid[~mask]
interp = interpolate.griddata((x1, y1), newarr.ravel(), (xx, yy), method=method)
if method == "linear":
interp = interp_nans(interp, "nearest")
return interp
a_title = "Argo Data"
s_title = "Seal Data"
comb_title = "Argo+Seal Data"
interp_title = "Interpolated Argo+Seal Data"
In [137]:
resolution = 2
argo_pt_grid = load_netcdf_grid(target_var="n_point", with_depth = False, resolution = resolution)
seal_pt_grid = load_netcdf_grid(path="seal", target_var="n_point", with_depth = False, resolution = resolution)
In [5]:
combined_pt_grid = np.nansum((argo_pt_grid, seal_pt_grid), axis=0)
combined_pt_grid[np.isnan(argo_pt_grid) & np.isnan(seal_pt_grid)] = np.nan
In [6]:
#plot_grid(argo_pt_grid, title="Number of vertical profiles in gridded argo dataset", cbtitle="Number of profiles", vmin=1,vmax=25, resolution = resolution)
plot_grid(argo_pt_grid, title=a_title + "\nNumber of vertical profiles", cbtitle="Number of profiles", vmin=1,vmax=25, resolution = resolution)
plot_grid(seal_pt_grid, title=s_title + "\nNumber of vertical profiles", cbtitle="Number of profiles", vmin=1,vmax=25, resolution = resolution)
plot_grid(combined_pt_grid, title=comb_title + "\nNumber of vertical profiles", cbtitle="Number of profiles", vmin=1,vmax=25, resolution = resolution)
In [13]:
resolution = 2
In [7]:
argo_temp_grid = load_netcdf_grid(resolution = resolution)
seal_temp_grid = load_netcdf_grid("seal", resolution = resolution)
In [8]:
combined_temp_grid = np.nanmean((argo_temp_grid, seal_temp_grid), axis=0)
In [9]:
title = a_title + "\n0-10 decibars mean temperature"
plot_grid(argo_temp_grid, title=title, resolution = resolution)
#title = "Seal data gridded at {}\n0-10 decibar mean temperature\nMean across all points = {}".format(1/resolution,sealmean)
title = s_title + "\n0-10 decibar mean temperature"
plot_grid(seal_temp_grid, title=title, resolution = resolution)
title = comb_title + "\n0-10 decibar mean temperature"
plot_grid(combined_temp_grid, title=title, resolution = resolution)
In [10]:
interp = interp_nans(combined_temp_grid)
In [11]:
#title = "Argo+Seal data gridded at {}\n0-10 decibar mean temperature\nMean across all points = {}".format(1/resolution,interpmean)
title = interp_title + "\n0-10 decibar mean temperature"
plot_grid(interp, resolution = resolution, title=title)
In [19]:
argo_max_depth = load_netcdf_grid(target_var = "depth", resolution = resolution)
seal_max_depth = load_netcdf_grid("seal", target_var = "depth", resolution = resolution)
In [20]:
plot_grid(argo_max_depth, resolution = resolution, title="Argo data Max depth\n", cbtitle="Depth (dbar)")
plot_grid(seal_max_depth, resolution = resolution, title="Seal data Max depth\n", cbtitle="Depth (dbar)", vmax=2000)
In [177]:
for month in range(1, 13):
try:
argo_temp = np.load(f"argo_temp_grid_{month}.npy")
seal_temp = np.load(f"seal_temp_grid_{month}.npy")
combined_temp = np.load(f"combined_temp_grid_{month}.npy")
except:
argo_temp = load_netcdf_grid(filter_month = month, with_depth=True)
seal_temp = load_netcdf_grid("seal", filter_month = month, with_depth=True)
combined_temp = np.nanmean((argo_temp, seal_temp), axis=0)
np.save("argo_temp_grid_{}".format(month), argo_temp)
np.save("seal_temp_grid_{}".format(month), seal_temp)
np.save("combined_temp_grid_{}".format(month), combined_temp)
interp_temp = interp_nans(combined_temp[:,:,0])
argomean = np.round(np.nanmean(argo_temp), 2)
sealmean = np.round(np.nanmean(seal_temp), 2)
combmean = np.round(np.nanmean(combined_temp), 2)
interpmean = np.round(np.nanmean(interp_temp), 2)
# #title = a_title + "\n0-10 decibar mean temperature for {}\nMean across all points = {}".format(calendar.month_abbr[month], argomean)
title = a_title + "\n0-10 decibars mean temperature for {}".format(calendar.month_abbr[month])
#plot_grid(argo_temp[:,:,0], title=title)
# #title = "Seal data\n0-10 decibar mean temperature for {}\nMean across all points = {}".format(calendar.month_abbr[month], sealmean)
title = s_title + "\n0-10 decibar mean temperature for {}".format(calendar.month_abbr[month])
#plot_grid(seal_temp[:,:,0], title=title)
title = comb_title + "\n0-10 decibar mean temperature for {}".format(calendar.month_abbr[month])
# title = "Argo+Seal data\n0-10 decibar mean temperature for {}\nMean across all points = {}".format(calendar.month_abbr[month], combmean)
#plot_grid(combined_temp[:,:,0], title=title)
title = interp_title + "\n0-10 decibar mean for {}".format(calendar.month_abbr[month])
#plot_grid(interp_temp, title=title)
In [179]:
resolution = 1
for season in ["Summer", "Winter"]:
try:
#argo_pts_grid = np.load(f"argo_pts_grid_{season}.npy")
argo_temp = np.load(f"argo_temp_grid_{season}.npy")
seal_temp = np.load(f"seal_temp_grid_{season}.npy")
combined_temp = np.load(f"combined_temp_grid_{season}.npy")
except:
argo_temp = load_netcdf_grid(filter_season = season, with_depth=True, resolution = resolution)
#added number of point for summer season
#argo_pts_grid = load_netcdf_grid(target_var="n_point", filter_season = season, with_depth = True, resolution = resolution)##
seal_temp = load_netcdf_grid("seal", filter_season = season, with_depth=True, resolution = resolution)
combined_temp = np.nanmean((argo_temp, seal_temp), axis=0)
#np.save("argo_pts_grid_{}".format(season), argo_pts_grid) ##
np.save("argo_temp_grid_{}".format(season), argo_temp)
np.save("seal_temp_grid_{}".format(season), seal_temp)
np.save("combined_temp_grid_{}".format(season), combined_temp)
interp_temp = interp_nans(combined_temp[:,:,0])
argomean = np.round(np.nanmean(argo_temp), 2)
sealmean = np.round(np.nanmean(seal_temp), 2)
combmean = np.round(np.nanmean(combined_temp), 2)
interpmean = np.round(np.nanmean(interp_temp), 2)
#title = "Argo data 0-10 decibar mean temperature for {}\nMean across all points = {}".format(season, argomean)
title = a_title + "\n0-10 decibar mean temperature for {}".format(season)
plot_grid(argo_temp[:,:,0], title=title, resolution = resolution)
title = s_title + "\n0-10 decibar mean temperature for {}".format(season)
plot_grid(seal_temp[:,:,0], title=title, resolution = resolution)
title = comb_title + "\n0-10 decibar mean temperature for {}".format(season)
plot_grid(combined_temp[:,:,0], title=title, resolution = resolution)
title = interp_title + "\n0-10 decibar mean for {}".format(season)
plot_grid(interp_temp, title=title, resolution = resolution)
#plot_grid(argo_pts_grid, title="Number of points in gridded argo dataset", cbtitle="Number of points", vmin=1,vmax=25, resolution = resolution)
#title="Argo data\nNumber of profiles for {}".format(season) ##
#plot_grid(argo_pts_grid, title=title, cbtitle="Number of profiles", vmin=1,vmax=25, resolution = resolution)##
In [31]:
resolution = 2
try:
argo_temp_grid_withdepth = np.load("argo_temp_grid_withdepth.npy")
seal_temp_grid_withdepth = np.load("seal_temp_grid_withdepth.npy")
combined_temp_grid_withdepth = np.load("combined_temp_grid_withdepth.npy")
except:
argo_temp_grid_withdepth = load_netcdf_grid(with_depth=True, resolution = resolution)
seal_temp_grid_withdepth = load_netcdf_grid("seal", with_depth=True, resolution = resolution)
combined_temp_grid_withdepth = np.nanmean((argo_temp_grid_withdepth, seal_temp_grid_withdepth), axis=0)
np.save("argo_temp_grid_withdepth", argo_temp_grid_withdepth)
np.save("seal_temp_grid_withdepth", seal_temp_grid_withdepth)
np.save("combined_temp_grid_withdepth", combined_temp_grid_withdepth)
In [145]:
dbars = [0, 5, 10, 15, 20, 30, 50, 75, 100, 125, 150, 175, 190, 199]
tmax = 6
tmin = -3
In [18]:
for db in dbars:
mtf = "mean temperature field"
#Argo data\n0-10 decibar mean temperature for {}\nMean across all points = {}".format(calendar.month_abbr[month], argomean)
dbs = "{}-{} dbar ".format(db*10, (db+1)*10)
# dbs = "between {}-{} dbar".format(db*10, (db+1)*10)
plot_grid(argo_temp_grid_withdepth[:,:,db], title=a_title + "\n" + dbs + mtf, resolution=resolution, vmin=tmin, vmax=tmax)
plot_grid(seal_temp_grid_withdepth[:,:,db], title=s_title + "\n" + dbs + mtf, resolution=resolution, vmin=tmin, vmax=tmax)
plot_grid(combined_temp_grid_withdepth[:,:,db], title=comb_title + "\n" + dbs + mtf, resolution=resolution, vmin=tmin, vmax=tmax)
interp = interp_nans(combined_temp_grid_withdepth[:,:,db])
plot_grid(interp, title=interp_title + "\n" + dbs + mtf, resolution=resolution, vmin=tmin, vmax=tmax)
In [180]:
fig = plt.figure(figsize=(15,15))
m = Basemap(projection='spstere',boundinglat=-55,lon_0=180,resolution='l')
#m.drawcoastlines()
#m.fillcontinents(color='black',lake_color='aqua')
# draw parallels and meridians.
m.drawparallels(np.arange(-60, 0, 20))
m.drawmeridians(np.arange(-180, 181, 20), labels=[1,1,0,1])
#m.drawmapboundary(fill_color='aqua')
m.etopo()
plt.title("Bathymetry south of -60S")
Out[180]:
In [28]:
def plot_cross_slope(grid, title="Temperature at 179° W", lon=-179, cbtitle="Temperature °C", resolution = 1):
fig, ax = plt.subplots(1, 1, figsize=(20,10))
y = np.arange(-60, -75, -1/resolution)
z = np.arange(0, 2000, 10)
levels = np.arange(np.nanmin(grid), np.nanmax(grid), .01)
im = ax.contourf(y, z, grid.T, cmap="jet", levels = levels, extend="both")
ax.set_xlabel("° Latitude")
ax.set_ylabel("Depth (dbar)")
ax.set_title("{} at {} resolution".format(title, 1 / resolution))
cb = fig.colorbar(im)
cb.set_label(cbtitle, rotation=270)
db = oceansdb.ETOPO()
y = np.arange(-60, -75, -.01)
h = -db['topography'].extract(lat=y, lon=lon)["height"]
#ax.plot(y, h)
ax.fill_between(y, 5000, h, color='black')
ax.set_ylim(2000, 0)
ax.set_xlim(-60, -74.5)
plt.show()
filtered = combined_temp_grid_withdepth[:,2,:]
interp = interp_nans(filtered)
plot_cross_slope(filtered, resolution = resolution)
plot_cross_slope(interp, title="Interpolated Temperature at 179° W", resolution = resolution)
In [29]:
try:
argo_sal_grid_withdepth = np.load("argo_sal_grid_withdepth.npy")
seal_sal_grid_withdepth = np.load("seal_sal_grid_withdepth.npy")
combined_sal_grid_withdepth = np.load("combined_sal_grid_withdepth.npy")
except:
argo_sal_grid_withdepth = load_netcdf_grid(with_depth=True, target_var="sal", resolution = resolution)
seal_sal_grid_withdepth = load_netcdf_grid("seal", with_depth=True, target_var="sal", resolution = resolution)
combined_sal_grid_withdepth = np.nanmean((argo_sal_grid_withdepth, seal_sal_grid_withdepth), axis=0)
np.save("argo_sal_grid_withdepth", argo_sal_grid_withdepth)
np.save("seal_sal_grid_withdepth", seal_sal_grid_withdepth)
np.save("combined_sal_grid_withdepth", combined_sal_grid_withdepth)
In [31]:
fig, ax = plt.subplots(1, 1, figsize=(20,10))
print(argo_sal_grid_withdepth.shape)
seal_counts = np.sum(~np.isnan(seal_sal_grid_withdepth), axis=2)
argo_counts = np.sum(~np.isnan(argo_sal_grid_withdepth), axis=2)
print(np.argwhere(seal_counts + argo_counts >= 300))
depth = np.arange(0, 2000, 10)
ax.set_xlabel("Salinity")
ax.set_ylabel("Depth (dbar)")
ax.set_title("Salinity profile")
ax.set_ylim(2000, 0)
ax.plot(argo_sal_grid_withdepth[11, 466, :], depth)
ax.plot(seal_sal_grid_withdepth[11, 466, :], depth)
ax.legend(["argo", "seal"])
Out[31]:
In [32]:
for season in ["Summer", "Winter"]:
argo_temp = np.load("argo_temp_grid_{}.npy".format(season))
seal_temp = np.load("seal_temp_grid_{}.npy".format(season))
fig, ax = plt.subplots(1, 1, figsize=(20,10))
depth = np.arange(0, 2000, 10)
ax.set_xlabel("{} temperature °C".format(season))
ax.set_ylabel("Depth (dbar)")
ax.set_title("{} temperature °C profile".format(season))
ax.set_ylim(2000, 0)
ax.plot(argo_temp[6, 233, :], depth)
ax.plot(seal_temp[6, 233, :], depth)
ax.legend(["argo", "seal"])
In [33]:
dbars = [0, 5, 10, 15, 20, 30, 50, 75, 100, 125, 150, 175, 190, 199]
vmax = 35
vmin = 33.5
for db in dbars:
msf = "mean salinity field"
dbs = "{}-{} dbar ".format(db*10, (db+1)*10)
plot_grid(argo_sal_grid_withdepth[:,:,db], title=a_title + "\n" + dbs + msf , resolution=resolution, vmin=vmin, vmax=vmax)
plot_grid(seal_sal_grid_withdepth[:,:,db], title=s_title + "\n" + dbs + msf, resolution=resolution, vmin=vmin, vmax=vmax)
plot_grid(combined_sal_grid_withdepth[:,:,db], title=comb_title + "\n" + dbs + msf, resolution=resolution, vmin=vmin, vmax=vmax)
interp = interp_nans(combined_sal_grid_withdepth[:,:,db])
plot_grid(interp, title=interp_title + "\n" + dbs + msf, cbtitle="Salinity", resolution=resolution, vmin=vmin, vmax=vmax)
In [34]:
def dens_smow(T):
'''
Function calculates density of standard mean ocean water (pure water) using EOS 1980.
INPUT: T = temperature [degree C (ITS-90)]
OUTPUT: dens = density [kg/m^3]
'''
a0 = 999.842594;a1 = 6.793952e-2;a2 = -9.095290e-3;a3 = 1.001685e-4;a4 = -1.120083e-6;a5 = 6.536332e-9
T68 = T * 1.00024
dens = (a0 + (a1 + (a2 + (a3 + (a4 + a5*T68)*T68)*T68)*T68)*T68)
return dens
def dens0(S, T):
'''
Function calculates density of sea water at atmospheric pressure
USAGE: dens0 = dens0(S,T)
DESCRIPTION:
Density of Sea Water at atmospheric pressure using
UNESCO 1983 (EOS 1980) polynomial.
INPUT: (all must have same dimensions)
S = salinity [psu (PSS-78)]
T = temperature [degree C (ITS-90)]
OUTPUT:
dens0 = density [kg/m^3] of salt water with properties S,T,
P=0 (0 db gauge pressure)
'''
assert S.shape == T.shape
T68 = T * 1.00024
# UNESCO 1983 eqn(13) p17.
b0 = 8.24493e-1;b1 = -4.0899e-3;b2 = 7.6438e-5;b3 = -8.2467e-7;b4 = 5.3875e-9
c0 = -5.72466e-3;c1 = +1.0227e-4;c2 = -1.6546e-6
d0 = 4.8314e-4
dens = (dens_smow(T) + (b0+ (b1+(b2+(b3+b4*T68)*T68)*T68)*T68)*S + (c0+(c1+c2*T68)*T68)*S*np.sqrt(S)+ d0*(S**2) )
return dens
In [36]:
#Cell added on 21/10/2019 - still testing
# adding function for TS Diagram
dens = dens0(argo_sal_grid_withdepth, argo_temp_grid_withdepth)
# print(combined_temp_grid_withdepth[0])
#shape of dens - (30,720,200) = (latitude, longitude, vertical layers)
cmap = plt.get_cmap('jet_r');cmap.set_bad(color='white')
fig = plt.figure(4,figsize=(15,15)) #TS-diagram of Argo profiles in the Ross Sea Slice
plt.scatter(argo_sal_grid_withdepth[0:10],argo_temp_grid_withdepth[0:10],s = 2,c=dens[0:10],cmap = cmap, edgecolors='none',alpha=0.8)
plt.xlim(33.8,35);plt.ylim(-2.5,5)
plt.clb = plt.colorbar();plt.clb.ax.set_title('Density (kg/m$^3$)')
plt.xlabel('Pratical salinity',fontsize = 18);
plt.ylabel('Temperature ($^oC$)',fontsize = 18)
#ax.set_xlabel('Pratical salinity',fontsize = 18);
#ax.set_ylabel('Temperature ($^oC$)',fontsize = 18)
plt.title('T-S diagram, Argo profiles',fontsize = 20)
# plt.axis('tight')
plt.show()
# for i in range(0, len(FILES)):
# read = Dataset(source + count_60[j][1] + '/profiles/' + FILES[i], mode = 'r')
# Lat = read.variables['LATITUDE'][:]
# Lon = read.variables['LONGITUDE'][:]
# if Lon[0] < -140, Lon[0] > 160:
# Temp = read.variables['TEMP'][:]
# Psal = read.variables['PSAL'][:]
# sigma = dens0(Psal[0],Temp[0])
# scatter(Psal[0],Temp[0],s = 2,c=sigma,cmap = cmap, edgecolors='none',alpha=0.8)
# xlim(33.8,35);ylim(-2.5,5)
# clb = colorbar()#;clb.ax.set_title('Density (kg/m$^3$)')
# axis('tight');grid('on',alpha=0.2)
# # savepath = os.path.join('C:/Users/fvie285/Desktop/PhD_Project/Data/figures/Argo_South_of_50/aoml_figures/TS_diagram_RSS.png')
# # savefig(savepath, bbox_inches='tight', dpi=300)
# show()
In [37]:
dens = dens0(combined_sal_grid_withdepth, combined_temp_grid_withdepth)
dbars = [30, 50, 100, 150, 199]
vmax = 1028
vmin = 1026.5
for db in dbars:
dbs = "from {}-{} dbar".format(db*10, (db+1)*10)
interp = interp_nans(dens[:,:,db])
print(np.nanmin(interp), np.nanmax(interp))
plot_grid(interp, title="Interpolated gridded seal+argo density data " + dbs, cbtitle="Density", resolution=resolution, vmin=vmin, vmax=vmax)
In [38]:
plot_cross_slope(dens[:,2,:], title="Density at 179° W", cbtitle="Density", resolution=resolution)
plot_cross_slope(interp_nans(dens[:,2,:]), title="Interpolated density at 179° W", cbtitle="Density", resolution=resolution)
In [189]:
resolution = 2
years = [2007, 2011, 2015]
try:
grids_per_year = np.load("grids_per_year.npy", allow_pickle=True)
except:
grids_per_year = {}
for y in years:
for s in ["Summer", "Winter"]:
start = datetime(y,1,1)
end = datetime(y + 4,1,1)
argo_temp_grid_withdepth_yearfiltered = load_netcdf_grid(with_depth=True, resolution = resolution, filter_start = start, filter_end = end, filter_season = s)
seal_temp_grid_withdepth_yearfiltered = load_netcdf_grid("seal", with_depth=True, resolution = resolution, filter_start = start, filter_end = end, filter_season = s)
combined_temp_grid_withdepth_yearfiltered = np.nanmean((argo_temp_grid_withdepth_yearfiltered, seal_temp_grid_withdepth_yearfiltered), axis=0)
if y not in grids_per_year:
grids_per_year[y] = {
"Winter": {},
"Summer": {}
}
grids_per_year[y][s] = {
"argo": argo_temp_grid_withdepth_yearfiltered,
"seal": seal_temp_grid_withdepth_yearfiltered,
"combined": combined_temp_grid_withdepth_yearfiltered
}
np.save("grids_per_year", grids_per_year)
In [168]:
plt.rcParams['figure.facecolor']='white'
db = 50
for i, y in enumerate(years):
for s in ["Summer", "Winter"]:
argo_temp_grid_withdepth_yearfiltered = grids_per_year[y][s]["argo"]
seal_temp_grid_withdepth_yearfiltered = grids_per_year[y][s]["seal"]
combined_temp_grid_withdepth_yearfiltered = grids_per_year[y][s]["combined"]
title = f"at {db * 10}dbar from the year {y}-{y + 4} in {s}"
plot_grid(argo_temp_grid_withdepth_yearfiltered[:,:,db], title=a_title + "\n" + title, resolution=resolution, vmin=tmin, vmax=tmax)
plot_grid(seal_temp_grid_withdepth_yearfiltered[:,:,db], title=s_title + "\n" + title, resolution=resolution, vmin=tmin, vmax=tmax)
plot_grid(combined_temp_grid_withdepth_yearfiltered[:,:,db], title=comb_title + "\n" + title, resolution=resolution, vmin=tmin, vmax=tmax)
interp = interp_nans(combined_temp_grid_withdepth_yearfiltered[:,:,db])
plot_grid(interp, title=interp_title + "\n" + title, resolution=resolution, vmin=tmin, vmax=tmax)
if i > 0:
delta = argo_temp_grid_withdepth_yearfiltered[:,:,db] - grids_per_year[years[i - 1]][s]["argo"][:,:,db]
print(np.nanmin(delta), np.nanmean(delta), np.nanmax(delta))
plot_grid(delta, title=comb_title + "\n" + f"Delta from {years[i - 1]} - {y} in {s} at {db}dbar", resolution=resolution)
In [171]:
delta = grids_per_year[2015]["Winter"]["argo"][:,:,db] - grids_per_year[2011]["Winter"]["argo"][:,:,db]
print(np.nanmax(delta))
print(np.where(delta == np.nanmax(delta)))
delta[10,636]
print(10/2, 636 /2 - 180)
In [172]:
np.argo["lat"]
Out[172]: